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In the framework of density functional theory a formalism to describe electronic transport in the 
steady state is proposed which uses the density on the junction and the steady current as basic 
variables. We prove that, in a finite window around zero bias, there is a one-to-one map between 
the basic variables and both local potential on as well as bias across the junction. The resulting 
Kohn-Sham system features two exchange-correlation (xc) potentials, a local xc potential and an 
xc contribution to the bias. For weakly coupled junctions the xc potentials exhibit steps in the 
density-current plane which are shown to be crucial to describe the Coulomb blockade diamonds. 

At small currents these steps emerge as the equilibrium xc discontinuity bifurcates. The formalism 
is applied to a model benzene junction, finding perfect agreement with the orthodox theory of 
Coulomb blockade. 

PACS numbers: 05.60.Gg, 31.15.ee, 71.15.Mb, 73.63.-b 


Engineering electrical transport through molecular 
junctions is a mandatory passage toward the miniatur¬ 
ization and speeding up of device components ffll- 
As a systematic experimental characterization of every 
synthesizable molecule is impractical, reliable theoretical 
methods are of utmost importance to progress. Density 
Functional Theory (DFT) has emerged as the method of 
choice due to the chemical complexity of the junctions 0 - 
0. Nevertheless, to date there exists still no DFT scheme 
to deal with the ubiquitous Coulomb Blockade (CB) phe¬ 
nomenon. CB stems from the interplay between quantum 
confinement and Coulomb repulsion, and it leaves clear 
fingerprints in the measured conductances. 

As recognized by several authors the disconti¬ 

nuity of the exchange-correlation (xc) potential plays a 
pivotal role in blocking the electron density at integer val¬ 
ues. Although at low temperatures and for single-channel 
junctions this key xc feature yields also accurate DFT 
conductances [131 - 118] . we showed that at finite temper¬ 
ature and/or bias the exact DFT conductance does not 
exhibit any CB signature m- In fact, static DFT misses 
the dynamical xc bias corrections [T0H22] predicted by 
Time-Dependent (TD) DFT [53], the proper framework 
within which to formulate a theory of quantum trans¬ 
port [19|, 20j. Only recently a dynamical xc correction 
has been proposed El. but its applicability is limited to 
the CB regime at zero bias. 

In this Letter we put forward a steady-state DFT, 
henceforth named i-DFT, whose xc potentials (gate and 
bias) are functionals of the molecular density and the 
steady-state current. The i-DFT framework is suited 
to study the finite-bias and finite-temperature conduc¬ 
tance as function of external gate and applied bias, and 
it generalizes standard DFT in equilibrium. Unlike Cur¬ 
rent DFT [23] (CDFT), i-DFT applies out of equilibrium 
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FIG. 1. Schematic illustration of a molecular junction at¬ 
tached to a left and right electrode. The system is subject to 
a bias voltage V (solid line) and a gate potential iy;(r). The 
grey area is the molecular region defined in the main text. 

and unlike TD(C)DFT [251125] the i-DFT functional has 
no memory. The empirical Landauer+DFT approach to 
transport [3HZ] is recovered as an approximation to i- 
DFT with the xc bias set to zero. As we shall see this 
approximation is too severe in the CB regime. 

Through a reverse-engineering procedure we show that 
the well-known discontinuities of the DFT xc potential at 
integer particle number N bifurcate as the current / starts 
flowing. The xc gate and xc bias exhibit an intricate and, 
at first sight, inexplicable pattern of intersecting steps in 
the N-I plane. We recognize, however, that every inter¬ 
section occurs at the plateau values of N and I in a CB 
diamond. This “duality” between intersections and CB 
plateaus in current and particle number is an exact prop¬ 
erty of the xc potentials of nonequilibrium open systems 
and the fundamental ingredient to extend DFT transport 
calculations to finite bias. 

i-DFT We consider a current-carrying molecular junc¬ 
tion attached to a left (L) and right ( R ) electrode in the 
steady state, see Fig. [T] In addition to the nuclear po¬ 
tential v n (r) the electrons are subject to an external bias 
vs(r) generated by a battery and to an external gate 




2 


vc(r) that vanishes deep inside the electrodes. The clas¬ 
sical potential i;„(r)+i;B(r)+ujj(r), vh being the Hartree 
potential, in the far L ( R ) region differs by a uniform shift 
V/2 (— V/2) from its equilibrium value, V being the ap¬ 
plied bias. m Since a change of the external gate does 
not affect the value of V we can unambiguously calculate 
the steady-state current and density by specifying v n (r), 
vg(t) and V. [T9 1 f20| 

Let us select an arbitrary finite region of space 7Z 
around the molecule, henceforth named molecular re¬ 
gion, and write the nuclear/gate potential as v n /c{ r) = 

C/gM+^/gM where v n/c( v ) = v n/c{r) for r <E 7^ and 
zero otherwise. Similarly, we write the electronic density 
as n( r) = n ln (r) + n out (r) where n ln (r) = n(r) for r £ 7Z 
and zero otherwise. We advance that for practical ap¬ 
plications a convenient choice of the molecular region is 
the one for which in the electrodes (outside 1Z) we have 
vg — 0 and v n + vh weakly dependent on the chemi¬ 
cal structure of the junction. In this case n out (r) is not 
affected by a change in vq, a condition used in all DFT- 
based quantum transport calculations. However, the for¬ 
mal results contained in this section are independent of 
the choice of 7 Z. 

The foundation of i-DFT rests on the one-to-one 
correspondence between the two pairs (u ln (r),V) and 
(n m (r),J). The first pair consists of the molecular po¬ 
tential u ln (r) = ujf (r) + Vq (r) and the bias V whereas 
the second pair consists of the molecular density ri ln (r) 
and the steady current I. 

Theorem- For any finite temperature and at fixed 
outer potential v out (r) = t>£ Ut (r) + u° ut (r) the map 
(v ln (r), V) —> (n ln ( r), I) is invertible in a finite bias win¬ 
dow around V = 0. 

Proof.- To prove the theorem we show that the Jaco¬ 
bian 


Jy= o — Det 


<5?z ln (r)/<5-G ln (r') dn ln {r)/dV 
6I/6v in { V) dl/dV 


v=o 


( 1 ) 


is nonvanishing (we are of course working under the 
physically reasonable assumption that n ln (r) and I are 
continuously differentiable, which also implies that Jy 
is continuous in V = 0). The block y ln (r,r') = 
<5n m (r)/Ju m (r , )|y = o is the static equilibrium density re¬ 
sponse function of the contacted system with r, r' in the 
molecular region, whereas G = dl/dV\y= 0 is the zero- 
bias conductance. We start by showing that y ln is in¬ 
vertible for any molecular region 7 Z. The equilibrium re¬ 
sponse function % ln can be calculated using leads of finite 
length L and then taking the limit L —>• oo. Let {|\Eq}} 
be a complete set of many-body eigenstates of the equi¬ 
librium contacted system with energy Ej and number of 
particles TV*. At temperature 1//3 and chemical potential 


/i the Lehmann representation of y ln reads [28] 

X“(r,r') = ^ 
ij 

(2) 

with Z the partition function, f lij = Ei — Ej the energy 
difference, rj a positive infinitesimal to set to zero after 
the limit L —» oo, and fij( r) = (^j|n(r)|^j) — Sijn(r) the 
excitation amplitudes. We define = j n drfij (r)t(r) 
where t( r) is a test function. Notice that fij has a well 
defined limit for L —> oo since the integral is over the 
finite domain 7 Z. Proving the invertibility of y ln is equiv¬ 
alent to proving that 

I[ drdr't(r)^(ry)t(r') = \ ^ 

x (e~^ - e-W) e^ ± 0 (3) 

for any test function t(r). To this end we observe that for 
every Ei Ej we have f \j 0 and (e~P Ei — e~ l3Ej ) ^ 0. 

Hence the left hand side of Eq. © is non-positive and 
can be zero only provided that Tij = 0 for every Ei 
Ej. However, this latter circumstance implies that the 
Hamiltonian and the operator T = drn(r)t(r) can be 
diagonalized simultaneously, an absurdum for any test 
function t( r). Therefore Eq. © holds true and y m is 
invertible. Similarly, from the Lehmann representation 
of the zero-bias conductance [29] 


/u( r )/ii( r ') 


n 2 . 
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one finds the intuitive result G > 0. In Eq. Q 7 tJ = 
(T, |T|j) with I the longitudinal current operator. The 
crucial observation to end the proof is that at zero bias a 
variation of u ln does not induce a steady current, hence 
<5//<5u in (r')|y =0 = 0. We conclude that J 0 = Det[,\ in ]G < 
0 for all v ln . Since Jy is a continuous function of V 
around V = 0, there exists a finite interval (depending 
on v ln ) around V = 0 for which Jy < 0. In this domain 
the map (u ln (r),T) —► (?r ln (r),J) is invertible. 

An interesting consequence of the i-DFT theorem is 
that at zero bias (and hence at zero current) it generalizes 
standard equilibrium DFT at finite temperatures |3U] to 
portions of an interacting system. In fact, the i-DFT the¬ 
orem implies that two potentials v and v' differing only 
in a region 7 Z generate two equilibrium densities n( r) and 
n'{ r) which are certainly different in 1Z (see also Ref. EU. 
It is worth observing that the zero-bias i-DFT does not 
suffer from the DFT problem of infinite systems [52] since 
the map involves only the density and the potential in JZ. 
Furthermore, since we are not interested in the density 
outside 1Z, no assumption on the analyticity of the den¬ 
sity in position space is needed [55] ■ Interestingly, the 
zero-bias i-DFT for infinite systems can easily be gen¬ 
eralized to the time domain too. Indeed, the boundary 
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term in the Runge-Gross proof [53] vanishes identically 
since v — v' is, by definition, zero outside 7Z. 

Henceforth we omit the superscript “in” in the molec¬ 
ular density and potential; thus n( r) and v(r) are always 
calculated in r G 1Z. Let (n, I) be the molecular density 
and current induced by the potentials (v, V) in an in¬ 
teracting junction. Under the usual assumption of non¬ 
interacting u-representability, the i-DFT theorem guar¬ 
antees that a pair of potentials (v a , V s ) which reproduces 
the same (n, I) in a noninteracting junction is unique. 
Following the Kohn-Sham (KS) procedure we define the 
xc bias and Hartree-xc (Hxc) gate as 

V xc [n, I] = V s [n,I ] - V[n,I ], (5a) 

^Hxc[«,d](r) = v s [n,I]{ r) - v[n,I}( r). (5b) 


The nomenclature Hxc gate instead of Hxc molecular po¬ 
tential is used for brevity to indicate that uhxc is nonzero 
only in the molecular region. The self-consistent KS 
equations then read (hereafter f = f 

n(r) = 2 [ f(u + s a V+ ^ c ) A a( r,w), (6a) 

a=L,R** 

7 = 2 E / f(eo + s a V + Vxc )s a T( u). (6b) 

a=L,R^ 


In Eqs. Q, f(u) = l/(e /3 i“ J_/i i + 1) is the Fermi func¬ 
tion whereas s_r/l = ±. We also defined the par¬ 
tial spectral function A a (r, u>) = (r|C/(a;)r Q (u;)<? 1 '(u;)|r), 
with Q and Tq, the KS Green’s function and broaden¬ 
ing matrices [53], and the transmission function T = 
Tr [C?(w)r L (o;)5 t (w)r /? (a;)]. In Eqs. ((g]) one should note 
the presence of the KS bias V s = V + V xc instead of 
the bare bias V; this comes from the mapping (n, I) O 
(v s ,V s ). Although the KS Eqs. Q for n and I are for¬ 
mally identical to the TDDFT expressions of Ref. [19| ISO, 
we stress that the i-DFT xc potentials depend on the 
steady-state molecular density and current while the 
TDDFT xc potential depend on the full history of the 
molecular and lead densities. The augmented local char¬ 
acter of the i-DFT xc potentials is in agreement with 
similar findings in TDCDFT [33 [35]. 

The simplifications brought about by i-DFT are es¬ 
pecially evident when considering the zero-bias conduc¬ 


tance. From Eq. (6b) we have 


G = G„ 


l + ^c 


dr 


SV XC dn{ r) 
8n(r) dV 


(7) 


v=o 


with G s = —2 J /'(w)T(w) the zero-bias KS conduc¬ 
tance. Since I = 0 is the solution of the KS equa¬ 
tions with V = 0 and since V xc [n, 0] = 0 (for other¬ 
wise there would be a current flowing in the system 
at zero bias, in contradiction to the theorem) we have 


5V xc [n, I\/8n{r)\ v =o = 5V xc [n, 0]/<5n(r) = 0. Therefore 


G = 


G s 


1 -G s 


dV x g 

dl 


1=0 


( 8 ) 


The i-DFT correction to G s is physically more trans¬ 
parent than the TDDFT correction involving the zero- 
momentum and zero-frequency limit of the xc kernel 

mmmm- 

Although i-DFT has been formulated in r-space the 
same theorem and KS procedure apply to tight-binding 
models by replacing r with a site or orbital index. To 
highlight the distinctive features of the i-DFT poten¬ 
tials in the CB regime, hence above the Kondo temper¬ 
ature, we apply the formalism to a junction described 
by the Constant Interaction Model (CIM) with Hamilto¬ 
nian H = + h U Y;icrjtja' where n ia 

is the occupation operator of the z-th level with spin 
a. [33 13H As the electron-electron interaction is confined 
to the molecular junction the Hartree potential vanishes 
in the leads. If we choose the region 1Z as the set of in¬ 
teracting levels then the i-DFT theorem states that there 
is a one-to-one map between the pair ({e,},U) and the 
pair ({nj,/). 

Anderson model The CIM with one level coupled 
to two featureless leads, r a (w) = y/2, is equivalent to 
the Anderson model. Due to spin-degeneracy the i-DFT 
potentials depend on the particle number N = n\^ + 
and current I. For £\ = v, these can be calculated from 


N = J[f(u>-V/2) + f{u> + V/2)]A{u>-v) (9a) 

I =\j [/(w - V/2) - f(to + V/2)]A( U - v), (9b) 

with A(oj) the interacting spectral function. Above the 
Kondo temperature TV a good approximation for the 
spectral function is A(u) = ^(w — U) + (1 — y)£(w) 
with £(lo) = 7 /(w 2 + 7 2 / 4 ) 139]. We verified (not shown) 
that N and I are in excellent agreement with the results 
of the Rate Equations [33 @D3 (RE), the orthodox theory 
of CB for weakly coupled systems. For the Anderson 
model the map (v, V) —» (N , I) is invertible for all V 
(infinite bias window) since the finite-bias Jacobian is 


Jv 


-2 7 


a_|_a_ 

2 — b + — 6_ 


( 10 ) 


where a± = f f(w±V/2)A(w — v) and b± = f f(ui)[£(w — 
v—U=FV/2)—£(u>—v=pV/2)]. As a± < 0 and — 1 < b± < 0 
we find J v < 0 HI] - 

To find the xc potentials in Eqs. ([5| we invert the map 
{v,V) —t ( N,I) of Eqs. ([9| for U V 0 and for U = 0. 
In terms of the variables w± = v ± V/2 the problem is 
separable since the map reads N^f2I/ r y = 2 f /(w) A (w — 
w±) = n(w±), which we solve by the bisection method. 
The values of n £ [0,2] and therefore the domain spanned 
by the particle number and current is |/| < ^N for N £ 
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FIG. 2. Hartree xc gate (top) and xc bias (bottom) of the 
Anderson model. Energies are in units of U. 


[0,1] and |/| < J(2 — N) for N £ [1,2]. The reverse 
engineered xc potentials i>hxc and V xc are shown in Fig. 
[2] and three observations arise, (i) The Hxc gate (xc bias) 
exhibits smeared steps of height U /2 (U) along the lines 
IV = lT I/'J. Interestingly, the DFT xc discontinuity of 
f Hxc [IV, 0] in N = 1 bifurcates as current starts flowing, 
(ii) The signs of V xc and I are opposite (in agreement with 
the results of Ref. [42]) and the derivative ^p|/ = o < 0, 
thus setting G s as the upper limit for the interacting 
zero-bias conductance [see Eq. ([8])]. In Ref. [TT]we found 
that the zero-bias TDDFT correction at IV = 1 can be 
expressed in terms of the Hxc gate as ~ 1/( 7 G S 
A comparison with Eq. ^ suggests the existence of a 
relation (at least in the Anderson model) between V xc 
and uhxc since for the two schemes to agree ^p|j = o = 
~~7 d d l N° |/=o ^ 1 /G s - (in) At finite bias \V\ < U and 
particle number N = 1 the CB prevents current flow in 
the interacting system. In the same situation the KS level 
pins to the chemical potential and I would be large if it 
were not for the counteraction of the xc bias V xc = — V. 
The Landauer+DFT approach [3HZ] would therefore fail 
dramatically as it misses xc bias corrections. 

We look for a parametrization of the xc potentials with 
the following properties: (a) at zero current V xc = 0 and 
uhxc reduces to the Hxc gate of Ref. [IT] (b) occurrence 
of smeared steps of height U /2 for uh xc (and U for I4c) 
at N = 1 7/7 and (c) the derivatives t ^f L \i = o and 

^r|j=o are related as discussed in observation (iii). The 
xc potentials 




^Hxc 


s=± 


1 -1— atan- 
7r 


N- 


- 1 


7 

W 


a N+iI-1 
V XC [N, I] = -U V - atan- 7 

Z —/ jy 

S=± 


w 


(11a) 

(lib) 


have the required properties. If we choose the fitting 


FIG. 3. Hartree xc gate (top) and xc bias (bottom) of the 
three-degenerate level CIM. The edges of two steps with pos¬ 
itive and negative slopes are highlighted with dashed lines. 
Energies are in units of U. 


parameter W ~ O.I 67 /U, i-DFT with the xc potentials of 
Eqs. produces self-consistent currents and densities 
which are almost indistinguishable from the interacting 
ones (not shown), including the density plateaus at | and 
| as well as the 2:1 ratio of the step heights in the current 
[151 . Remarkably, all this is achieved without breaking 
the spin symmetry. 

Notice also that in Ref. [ml there was no xc bias and 
the current did not reach a steady-state value since the 
fitting parameter W was set to zero (in this case no self- 
consistent solution of the steady-state equation exists). 

Multiple-level degenerate CIM Let us generalize 
the analysis to a CIM Hamiltonian with A4 degenerate 
levels £i = v coupled to featureless leads, T a ^j(u) = 
Sij^//2. Due to degeneracy urxc is uniform and depends 
only on N = n »<t and I. Above Tk the particle 

number and current of the interacting and noninteracting 
CIM can be calculated from the RE, and then inverted 
by an adaptation of the iterative bisection algorithm of 
Ref. 321 Finally, the i-DFT potentials can be obtained 
by subtraction as in Eq. ([ 5 ]). The map (y,V) —> (IV,/) 
is invertible for all V and the codomain is |/| < 7 N for 
N £ [0, M] and \I\ < 7 (2A4-IV) for N £ [M,2M]. The 
xc potentials are shown in Fig. [5] for A4 = 3. Like in the 
Anderson model the Hxc gate (xc bias) exhibits smeared 
steps of height U/2 ( U ) but the pattern of their edges 
is more complex. The equilibrium xc discontinuities 
of uhxc [IV, 0] at integer N bifurcate with IV-dependent 
slopes, and at every high-current intersection the slopes 
of the step edges change. 

In the attempt of disentangling the intricate pattern of 
discontinuity-lines we realized the existence of a duality 
between intersections in the xc potentials and plateaus 
in the particle number and current. From the RE of 
the degenerate CIM a plateau is uniquely identified by a 
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couple of integers (to, n ) with m,n = 0,..., 2A4; hence 
the number of plateaus is correctly given by ( 2M + l) 2 . 
In the (to, n)-plateau with n > to the probabilities P{q), 
q = to, ..., n of finding q particles are all identical and 

given by P-\q) = P~> m = £" =ro 

other probabilities vanish. The corresponding particle 
number and current are therefore 


whereas all 


N = N n > m = P n > m 3 


J—m 

n —1 


i = i, 


n>m — 


;p„> m M-n 


(12a) 


(12b) 


and 


Similarly, for n < m one finds TV = TV„< m = TV m >„ 

I = I n <m = —Im>n- The intersections in the xc poten¬ 
tials of Fig.[3]occur precisely at the points (TV n > m , / n > m ) 
and (TV„< TO , I n < m ). Knowledge of these points allows us 
to generalize the parametrization in Eqs. <0 


v (M) 

y Hxc 


U 


2M — 1 


M = T e e 


K=1 s=± L 
2M-1 


1 H— at an 

7r 


2 , A $(N,I) 


vM[N,I] = -U Y E“ atan 


K =1 s=± 


A 


w 

w 


(13a) 


(13b) 


where A^(iV, I) is the piece-wise linear function of TV 
and / which vanishes along the step edge passing through 
(AT, 0) and having positive (s = +) or negative (s = —) 
slopes (for examples see dashed lines in top panel of 
Fig. [3h. Here W is the same fitting parameter as in 
Eq.(|ll|). We verified (not shown) that the self-consistent 
solution of the KS equations with the xc-potentials of 


Eqs. (13) are in excellent agreement with the RE results. 


Using the analytic parametrization of the xc bias in 


Eq. (13b I we can calculate the zero-bias conductance 
from Eq. (|8|. The result is 


G_ 

g: 


1 


1 I 2UG S *r^2M — 1 2M-K + 1 ~t A' + l 

' 77rW 2-^jK —1 x+( N ^ K ) 2 


(14) 


The correction to G s is large for integer TV’s and Eq. (14) 


(_L_ 

\2M-1S 


2 UG S 

"ynW ^ 2J\A—N+1 


can be approximated as ~ 1/[1 + 
jyT_)]. Thus, the height of the conductance peaks de¬ 
pends on the number of electrons in the junction; the 
closer we get to half-filling the larger the height is. 

Multiple-level CIM and finite bias conductance 


The potentials in Eq. (13) are not suited to study junc¬ 


tions with nondegenerate levels as the dependence on 
the local occupations cannot be reduced to a depen¬ 
dence on TV only. New interesting aspects arise which 
are best illustrated in a HOMO-LUMO CIM with en¬ 
ergies €i = tn/L + v of degeneracy Mh/l■ Let N H / L 
be the particle number on the HOMO/LUMO level, 
and vh xc [N h ,N l ,I](H/L) and V XC [N H ,N L ,I] be the 


HOMO/LUMO Hxc gate and xc bias respectively. For 
Nl = 0 (empty LUMO) we have uhxc(-H) = vn xc (L) 
(uniform Hxc gate) and the i-DFT potentials are given 


by Eqs. (13) with M = Mh- Similarly for N H = 2Mh 


(full HOMO) the Hxc gate is uniform and the i-DFT po¬ 
tentials are again given by Eqs. (131 but with M = Ml- 
At zero current TV# ~ TV and Nl — 0 or TV# ~ 2 Mh 
and Nl — TV — 2A4# are the only physically realizable 
occupations. Hence, at least for small currents, uhxc is 
uniform [45] and can be parametrized by combining the 
Hxc gate of two CIMs with degeneracy M = Mh/l- 
This argument remains valid for an arbitrary number of 
levels; below we therefore show how to construct the i- 
DFT potentials for the general case. 

Let n = {ni.. .um} be the occupations of levels 
{1...A4}, A i P [n] the degeneracy of the p-tli largest 
occupation and T>[n] the number of distinct densities. 
For instance if A4 = 5 and n = {| ! 5 ; 2 ’ I ’ 3 1 t ^ ien 
Mi = 2, M 2 = 3 and V = 2. We further define 
A//[n] = 2 Y^q=i A4 9 [n] as the maximum number of par¬ 
ticles in the first (p— 1) degenerate levels (AT\ = 0). The 
degeneracies M p are used to construct the following i- 
DFT potentials 

VHxc|n, I] = Y W H^cf N) [ N ~ KM, 1 ] 

P= 1 

X>[n]-1 

+7 £ E 

p=l s=± 

T>[n\ 

v x cM i] = Y V ^ p[n]) l N ~ KM, 1] 

p =1 

CI " H s . TV+^J-A4 +1 [n] 


2 TV+ ^I-J\r p+1 [n} 

1 + - atan- - ——- 

7T W 


(15a) 


—U ^ ^ — atan- 


p=l s=d= 


W 


.(15b) 


where, again, W is defined as before. 

The dependence on the local occupations enters ex¬ 
clusively through the M p . At the joining points (TV = 
TV p +i[n] and 1 = 0) between two consecutive we 

add a discontinuity with slopes ± 2 / 7 . In fact, the slope 
of the lines delimiting the domain of the i-DFT potentials 
of a A4-fold degenerate CIM are independent of M, see 
Figs. [ 2 ] and [3j The i-DFT potentials in Eqs. [ fUj] ) reduce 
to those in Eqs. (13) for equal occupations rii = N/M 
and to the equilibrium DFT potentials of Ref. |TT| for 
1 = 0. Furthermore, they can easily be generalized 
to CIM with local interactions \ Uijfii a njcr' and 

level-dependent broadenings T a i j = Sij 7 ;/ 2 - 

To demonstrate the improvement of i-DFT over Lan- 
dauer+DFT we calculate the finite-bias differential con¬ 
ductance dl/dV of a benzene junction and benchmark 
the results against the RE. The benzene is described by 
a six-level CIM with U = 0.5 eV and e* = e? + v where 
e° = - £ ° = 5 .O 8 eV, = -e^ = -e° = -2.54 
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FIG. 4. Differential conductance (in units of the quantum of 
conductance) calculated from Landauer+DFT (top), i-DFT 
(middle) and RE (bottom), for a six-level CIM modelling ben¬ 
zene. Dashed red lines delimit the low-bias region where i- 
DFT and RE agree. 


eV. |46] In Fig. [4] we report the dl/dV in the three ap¬ 
proaches. Above Tk no Kondo plateau is expected. The 
Landauer+DFT scheme (top panel) does instead pro¬ 
duce Kondo plateaus at zero [TOUTS] and nonzero biases. 
Furthermore, it misses several dl/dV lines as compared 
to RE (bottom panel), thus providing a completely er¬ 
roneous description of CB. The i-DFT scheme (middle 
panel) correctly suppresses the spurious Landauer+DFT 
Kondo plateaus, and reproduces the RE trend of the CB 
peak-heights, see Eq. (14) and discussion below. At small 
bias, i-DFT captures all the dl/dV lines present in the 
RE approach while at higher bias some lines are missing. 
This latter fact is not surprising as our model xc poten¬ 
tials is designed to be accurate only at small currents. 


In general, the RE approach requires the solution of 
a linear system of size 4 M . In contrast, i-DFT requires 
the solution of A4 coupled nonlinear equations. If one 
aims only at describing the low-bias features correctly, in 
i-DFT the solution of two coupled equations is sufficient, 
independent of the size of the system. 

Conclusions We propose the i-DFT framework to 
calculate the steady density and current of interacting 
junctions at finite bias. i-DFT is based on the invert- 
ibility of the map between (n, I) and (v, V) in a finite 
bias window. Unlike the Landauer+DFT approach, i- 
DFT naturally leads to xc bias corrections, and in con¬ 
trast to TD(C)DFT, the i-DFT xc potentials are history- 
independent. We unveil the complex structure of the i- 
DFT xc potentials in the CB regime by reverse engineer¬ 
ing the occupations and current of a CIM as obtained 
from the RE. We found that the bifurcations of the xc 
discontinuity as current starts flowing are pivotal for the 
correct description of CB. Similarly to the equilibrium 
xc discontinuity of standard DFT the bifurcations are an 


intrinsic property of the i-DFT potentials, and are ex¬ 
pected to occur in the r-space formulation as well. 

We also find an efficient parametrization of the i-DFT 
potentials for small currents and use it to calculate the fi¬ 
nite bias conductance of a model benzene junction. Com¬ 
pared to Landauer+DFT, i-DFT shows a clear improve¬ 
ment, being able to reproduce all the (small bias) dl/dV 
lines of the RE approach. 
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